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Introduction 


Described  herein  are  the  key  findings  and  results  associated  with  the  project  entitled 
“Investigation  of  ELF  Signals  Associated  with  Mine  Warfare,  A  University  of  Idaho 
and  Acoustic  Research  Detachment  Collaboration,  Phase  One.”  Extremely  low 
frequency  (ELF)  electromagnetic  signals  are  used  by  enemy  combatants  to  detect  and, 
subsequently,  to  incapacitate,  by  means  of  surface  and  subsurface  mines,  naval  vessels. 
This  program  is  of  high  importance  to  the  Navy  -  particularly  since  ELF  signals  are 
one  of  the  primary  signature  emissions  of  the  Navy’s  proposed  electric  ship  fleet. 

In  principle,  the  questions  that  are  being  asked  in  this  investigation  are:  1)  once  an 
ELF  signal  is  generated,  how  far  does  it  propagate  and  still  be  detectable  and  2)  how 
can  such  signals  be  modeled,  excited  and  measured?  To  this  end,  the  scenario 
considered  is  one  in  which  an  ELF  source  of  the  electric  or  magnetic  kind  is  located  in 
or  above  water,  such  as  a  lake  or  ocean.  This  source  stimulates  an  ELF  signal  that  is 
free  to  propagate  in  the  water  and  air,  and  is  reflected  by  various  material  interfaces, 
say  between  the  water  and  air,  or  between  the  water  and  the  floor.  For  purposes  of 
experimental  demonstration,  the  investigation  focuses  on  the  scenario  of  ELF  sources 
and  signals  in  the  context  of  Lake  Pend  Oreille,  where  the  Acoustic  Research 
Detachment  (ARD,  Bayview,  Idaho)  is  located  and  entrusted  with  the  necessary  assets 
to  perform  validation  measurements. 

The  research  program  was  designed  with  two  major  thrusts:  Modeling  and 
experimentation.  The  modeling  thrust  was  coordinated  and  executed  by  the  University 
of  Idaho  (Ul),  Moscow,  Idaho;  the  experimentation  thrust  was  coordinated  and 
executed  by  ARD.  This  report  focuses  primarily  on  the  modeling  thrust.  A  separate 
report  from  ARD  will  be  issued  that  addresses  the  experimentation  thrust. 

ELF  Modeling 

Modeling  of  ELF  electromagnetic  signals  in  water  environments  can  be  accomplished 
either  by  means  of  direct,  analytical  solution  of  Maxwell’s  equations  or  by  numerical 
solutions  of  the  same.  The  former  is  attractive  for  purposes  of  gaining  insights  into  the 
physical  mechanisms  that  hinder  or  aid  the  propagation  of  ELF  signals.  The 
disadvantage  is  found  in  the  number  of  simplifying  assumptions  that  are  made  to 
bring  about  a  closed-form  solution.  A  numerical  solution  has  no  such  simplifying 
assumptions,  but  does  suffer  from  discretization  errors.  In  principle,  it  can  model  all 
of  the  physical  and  geometrical  features  of  the  domain  of  consideration.  The  price 
paid  for  doing  so,  however,  is  the  required  time  and  the  CPU/memory  resources 
needed  to  accomplish  the  task.  Data  visualization  and  management  are  other  issues 
that  need  to  be  addressed  when  working  with  large  data  sets  produced  by  numerical 
solvers.  The  positive  and  negative  tradeoffs  between  these  two  approaches  (i.e. 
analytical  vs.  numerical)  suggest  that  no  one  method  is  superior.  For  that  reason  the 
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U1  team  adopted  a  diverse  strategy  that  encompasses  many  different  approaches  in 
order  to  assure  a  positive  outcome  and  to  provide  deliverable  modeling 
methodologies. 

The  four  principle  techniques  that  were  considered  during  the  Phase  One  effort  were 
the  a)  Sommerfeld  Half-Space  (SHS)  method,  b)  Finite-difference,  time-domain 
method  (FDTD),  c)  High  Frequency  Structural  Simulator  (HFSS),  Finite-element 
method,  and  d)  the  quasi-static  method  (QS).  A  summary  of  these  methods  is  provided 
next.  Detailed  information  on  each  of  them  is  provided  in  the  attached  appendices. 

The  SHS  method  is  an  analytical  approach  that  assumes  that  all  interfaces  (say 
between  water  and  air,  or  between  water  and  floor)  are  planar  and  infinitely  extended. 
This  is  assumption  is  reasonably  valid  for  the  water-air  interface,  particularly  in  open 
water  regions  where  the  source  is  located  near  the  surface.  For  the  littoral  zones,  the 
method  may  fail,  particularly  when  electric  sources  are  used  to  excite  the  ELF  signals. 
By  assuming  that  the  interfaces  are  flat,  a  closed-form  solution  can  be  devised  that  is 
cast  in  terms  of  Fourier-Bessel  integrals.  These  integrals  can  be  evaluated  numerically 
and  rapidly  in  a  matter  of  seconds  on  any  desktop  machine.  Even  with  the  potential 
deficiency  of  treating  all  interfaces  as  planar,  the  SHS  method  is  attractive  as  a 
validation  tool  for  the  other  numerical  modeling  approaches.  For  example,  the  team 
used  the  SHS  method  to  validate  the  data  produced  by  the  FDTD  or  HFSS  methods 
(described  next)  when  these  numerical  methods  consider  the  same  layered  media 
problem  statement.  Professor  Robert  Olsen  of  Washington  State  University  (WSU)  is 
the  lead  investigator  of  the  SHS  method. 

The  FDTD  method  is  a  numerical  approach  that  discretizes  Maxwell’s  equations  in 
their  fundamental  form  using  a  staggered  grid  and  leap-frog  integrator.  This  method 
has  been  fully  vetted  in  the  open  literature  and  has  been  established  as  a  robust  way  of 
obtaining  accurate  simulation  data.  In  principle,  the  FDTD  method  accounts  for  all 
material  interfaces  and  material  inhomogeneities  by  assigning  permittivity, 
permeability  and  conductivity  values  along  edges  of  the  grid  elements.  Curvilinear 
boundaries  are  approximated  by  straight  line,  stair-stepped  boundaries.  For 
geometrical  features  that  are  significantly  less  than  a  wavelength,  such  stair-stepping 
causes  no  appreciable  errors  in  the  computed  data.  The  domain  of  interest  at  Lake 
Pend  Oreille  is  about  8  km  on  a  side;  the  lake  floor  at  its  deepest  point  is  about  335 
km.  Assuming  an  operating  frequency  of  100  Hz  and  a  water  conductivity  of  0.018 
S/m,  the  corresponding  skin  depth  is  375  m  and  the  wavelength  is  2.356  km;  for  air, 
the  wavelength  is  3,000  km.  Thus  the  domain  spans  a  fraction  of  a  wavelength  in  air 
but  about  3.4  wavelengths  (or  21.3  skin  depths)  in  water.  The  significant  disparity 
between  these  two  relative  sizes  potentially  introduces  computational  complexities. 
One  area  of  concern  that  was  addressed  in  this  project  was  the  proper  design  of  an 
absorbing  boundary  condition  (ABC)  or  perfectly  matched  layer  (PML)  that  will 
allow  an  open  physical  domain  to  be  truncated  into  a  finite  computational  domain. 
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Placement  of  this  ABC/PML  in  terms  of  wavelengths  is  critical  if  non-spurious 
reflections  are  to  be  avoided.  Typically,  the  farther  an  ABC  is  placed  away  from  a 
scattering  object,  the  better  it  performs.  However,  this  increases  the  computational 
domain  size  and  hence,  increases  the  simulation  time  and  CPU/memory  requirements 
of  the  computer.  PML’s  are  more  complex  to  implement,  but  can  be  placed  closer  to 
the  scattering  object  while  sustaining  excellent  absorption  characteristics.  Professor 
Dennis  Sullivan  of  the  University  of  Idaho  is  the  lead  engineer  of  the  FDTD  and  PML 
effort. 

HFSS  is  a  commercially  available  electromagnetic,  finite-element,  frequency-domain, 
numerical  solver  that  has  been  designed  by  Ansoft/Ansys  for  antenna  and  microwave 
circuit  applications.  One  question  that  was  asked  in  this  investigation  was  whether 
such  a  tool  could  be  used  to  predict  the  electromagnetic  propagation  characteristics  of 
an  ELF  signal  in  a  highly  conductive  environment.  We  believe  that  the  answer  is  a 
qualified  yes.  However,  being  that  the  solver  is  commercial  and  not  custom,  the  user 
does  not  have  direct  control  over  some  of  the  solver’s  functionality.  This  is  readily 
apparent  in  the  method  by  which  the  ELF  sources  are  constructed  in  the  solver's 
graphical  user  interface  (GUI).  Second,  some  of  its  interpolation  functions  cause  some 
errors  in  the  outputted  data.  Third,  learning  how  to  import  the  geometrical  features  of 
the  Lake  into  the  solver  and  optimizing  the  corresponding  computational  mesh  is 
needed.  Hence,  the  team  spent  considerable  amount  of  time  wrestling  with  these 
issues,  but  the  results  generated  to  date  seem  favorable.  Professor  Jeffrey  Young  of  the 
University  of  Idaho  is  the  lead  engineer  of  the  HFSS  effort. 

The  quasi-static  method  was  also  considered  given  the  ELF  nature  of  the 
electromagnetic  signals.  By  definition,  the  quasi-static  method  does  not  consider  any 
wavelike  mechanisms  in  Maxwell’s  equations;  it  assumes  that  the  field  lines  are  the 
same  as  the  static  field,  but  oscillating.  This  is  accomplished  by  neglecting  a)  electric 
displacement  currents  for  magnetic  sources  and  b)  magnetic  displacement  currents  for 
electric  sources.  By  doing  so,  simple  solutions  can  be  constructed  that  correlate  well 
with  other  more  advanced  solutions,  like  HFSS  and  FDTD.  Professor  Jeffrey  Young 
of  the  University  of  Idaho  is  the  lead  engineer  of  the  quasi-static  modeling  effort. 

Each  of  the  aforementioned  numerical  methods  requires  a  precise  understanding  of 
the  electrical  and  geometrical  features  of  the  Lake.  To  this  end,  the  team  spent 
considerable  amount  of  time  digitizing  graphical  topological  and  bathymetry  data 
over  an  8  km  by  8  km  region  of  the  Lake  that  is  centered  on  the  electromagnetic  array 
(EMA),  as  shown  in  Figure  1 . 
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Figure  1:  Electromagnetic  domain  of  interest. 

Renditions  of  the  cross-sectional  portion  of  the  Lake,  as  obtained  from  this 
discretization  process  and  starting  from  the  north  end  of  the  Lake,  are  shown  in 
Figures  2  through  5. 


Figure  2:  Cross-sectional  rendering  of  the  Lake;  far-north  end. 


Figure  3;  Cross-sectional  rendering  of  the  Lake;  middle-north  end. 


Figure  4:  Cross-sectional  rendering  of  the  Lake;  middle-south  end. 


Figure  5:  Cross-sectional  rendering  of  the  Lake;  far-south  end. 

For  these  cross-sections,  we  see  that  the  Lake  spans  about  6.5  km  and  has  a  depth  of 
about  375  m;  the  highest  land  feature  is  about  425  m  above  the  Lake.  Once  digitized 
and  stored  electronically,  each  of  the  previously  discussed  solvers  could  use  that  data 
in  their  computations. 

In  addition  to  precise  geometrical  data,  the  various  solvers  also  require  precise 
knowledge  of  the  conductivity  of  the  Lake  and  the  mud  at  the  bottom  of  the  Lake.  The 
U1  team  used  a  value  of  0.0 18  S/m  for  the  water  and  0.012  S/m  for  the  mud  floor.  The 
former  number  was  previously  measured  by  ARD;  the  latter  number  was  measured 
during  the  course  of  Phase  One  by  filling  a  PVC  tube  with  the  mud,  placing  electrodes 
on  the  ends  of  the  tube  as  caps  and  measuring  the  total  resistance  between  the 
electrodes.  By  knowing  the  geometrical  dimensions  of  the  tube,  we  were  able  to 
determine  the  conductivity  from  the  resistance  value.  As  for  the  value  of  the  dielectric 
permittivity  of  the  Lake,  this  was  not  deemed  essential,  since  displacement  currents  in 
the  Lake  are  virtually  insignificant  relative  to  the  conduction  currents. 

Although  detailed  information  pertaining  to  the  experimental  portion  of  this  project 
will  be  documented  by  ARD,  a  few  words  about  it  are  needed  to  set  the  context  of 
how  the  models  are  validated.  Most  of  the  validation  data  that  has  been  used  to  date 
was  gathered  in  December  2008  by  ARD.  In  one  set  of  experiments,  an  electric  source 
was  devised  by  mounting  two  metallic  plates  on  the  hull  of  a  boat  (see  Figure  6);  in 
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another,  one  plate  was  attached  to  the  boat,  which  also  towed  a  skiff  with  the  other 
metallic  plate  attached  to  it.  The  distance  between  the  boat  and  the  skiff  was  40 
meters. 


Figure  6:  Electrodes  on  the  source  boat. 

In  both  cases,  the  plates  were  connected  to  cables,  which  were  also  connected  to  a 
low-frequency,  high-powered  ELF  signal  source.  For  the  case  of  the  skiff',  an  ac 
electric  current  of  1.25  A  passed  through  the  plates  allowing  conduction  currents  to 
flow  in  the  water.  These  currents  stimulated  both  and  magnetic  fields  within  the  Lake. 
As  the  boat  motored  across  the  Lake’s  surface,  precise  GPS  data  was  logged  to 
determine  its  position  and  heading.  Over  five  days,  about  225  runs  were  executed 
consisting  of  85  noise/calibration  runs  and  140  experimental  runs.  The  experimental 
runs  for  phase  II  are  shown  in  Figure  7.  Many  of  the  boat  runs  were  redundant. 


Figure  7:  Experimental  boats  runs  conducted  by  ARD  in  December  2008. 


7 


The  grid  of  Figure  7  was  rotated  to  match  the  alignment  of  the  electromagnetic  array 
(EMA,  i.e.  field  sensing  unit)  thus  making  the  x-axis  shifted  10  degrees  positive  from 
that  of  easting.  The  alignment  of  the  sensors  can  randomly  shift  up  to  5  degrees 
during  the  measurement  process,  due  to  water  motion.  The  range  of  the  experiment 
was  about  6.7  km  in  the  x-direction  and  5.6  km  in  the  y-direction,  or  37.5  km  with 
the  EMA  as  the  origin.  Most  of  the  experiment  runs  were  perpendicular  or  transverse 
to  an  axis  of  the  EMA  and  some  were  at  arbitrary  angles  to  it. 

Of  these  140  runs,  about  25  of  them  have  been  processed.  The  average  boat  speeds  of 
these  runs  were  computed  from  the  GPS  and  time-stamping  data;  several  boat  speeds 
are  shown  in  Figure  8.  From  this  plot,  we  see  that  the  speeds  are  typically  between 
one  to  three  meters  per  second. 
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Figure  8:  Boat  speeds  for  various  runs. 


To  detect  the  fields,  a  series  of  electric  and  magnetic  sensors  were  connected  to  a  truss 
unit  called  the  electromagnetic  array  (EMA)  (see  Figures  9  and  10)  that  could  be 
raised  (or  lowered)  from  the  Lake  floor  to  the  surface  by  ARD’s  Target  Modeling 
Hauler  (TMH).  The  sensors  were  connected  to  a  data  acquisition  system  and  field 
strength  amplitudes  were  sampled  6,000  times  per  second,  resulting  in  a  Nyquist 
frequency  of  3  kHz. 
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Figure  9:  Electric  (.left)  and  magnetic  (right)  sensors. 


Figure  10:  Electromagnetic  Array 

Once  the  raw  experimental  data  was  gathered  and  recorded,  a  post-processing 
methodology  had  to  be  devised  that  could  extract  the  ELF  spectral  component 
associated  with  the  ELF  source.  This  processing  step  required  the  use  of  a  rectangular 
window  function  and  Fast  Fourier  Transform  (FFT)  techniques  operating  on  a  one 
second  sample  of  time-domain  data  (i.e.  6,000  samples).  The  other  part  of  the 
post-processing  step  required  the  manipulation  of  the  GPS  data  to  match  the 
coordinate  systems  of  the  FDTD,  UFSS  and  QS  models  with  the  experiment.  This 
step  is  not  necessarily  a  complicated  one,  but  care  does  need  to  be  exercised  to  insure 
that  the  various  coordinate  axes  are  parallel  and  share  a  common  origin.  Several 
MatLab  scripts  have  been  written  to  handle  both  the  FFT  and  GPS  processing  steps.  A 
typical  FFT  spectrum  from  the  experimental  data  is  shown  in  Figure  11.  The  100  Hz 
signal  frequency  is  clearly  seen,  but  is  surrounded  by  a  60  Hz  spectra  and  its  plethora 
of  harmonics.  The  noise  floor  is  seen  to  be  roughly  10  nV/m. 
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Figure  1 1:  Typical  spectrum  of  the  experimental  data. 

It  should  be  noted  that  the  experiment  was  conducted  by  moving  the  source  about  the 
lake  for  a  fixed  observation  point.  The  simulations  were  conducted  under  the 
assumption  that  the  source  is  fixed  and  the  observation  points  are  free  to  move.  Since 
the  relative  distance  between  the  source  and  the  observation  points  is  the  critical 
dimension  and  since  the  medium  is  reciprocal,  reciprocity  can  be  invoked  to  assure 
that  both  approaches  generate  identical  data. 

Key  Results 

In  the  attached  appendices,  several  specific  results  are  given  for  each  method. 
However,  many  of  these  methods  have  been  validated  by  each  other  and  by 
experimental  data.  These  key  results  are  presented  next. 

First  consider  Figure  12,  which  shows  the  six  test  validation  scenarios  and  boat  paths. 
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Figure  12:  Six  validation  scenarios  and  the  corresponding  boat  paths. 

Data  for  these  six  scenarios  are  shown  in  Figures  13  through  19.  A  summary  of  this 
data  is  provided  next. 

•  Figure  13:  This  corresponds  to  runs  4932,  4931  and  4926  for  a  40  m,  1.25  A, 
1,000  Hz  electric  source.  Data  from  the  quasi-static,  Sommerfeld  (i.e.  WSU), 
FDTD  and  HFSS  methods  for  the  x-component  of  the  electric  field  arc  plotted 
as  a  function  of  radial  distance  when  the  EMA  is  located  at  a  depth  of  8.23  m. 
All  but  the  quasi-static  method  align  very  well  with  the  experimental  results. 
However,  since  this  data  corresponds  to  distances  between  500  and  3,000 
meters,  which  span  multiple  wavelengths  at  1,000  Hz  (i.e.  A.=118  m),  the 
quasi-static  method  is  not  appropriate  for  this  test  case.  Even  so,  it  still 
predicts  the  correct  trend  and  is  only  off  by  a  factor  of  two. 

•  Figure  14:  This  corresponds  to  runs  3046,  3044  and  3042  for  a  40  m,  1.2  A, 
100  Hz  electric  source.  Data  from  the  quasi-static,  Sommerfeld  (i.e.  WSU), 
FDTD  and  HFSS  methods  for  the  x-component  of  the  electric  field  are  plotted 
as  a  function  of  radial  distance  when  the  EMA  is  located  at  a  depth  of  8.23  m. 
In  this  case,  all  methods  fail  to  give  the  correct  experimental  value,  but 
nevertheless  come  close  to  it.  Since  the  WSU  and  FDTD  methods  predict  the 
same  data,  we  surmise  that  the  error  is  associated  with  some  imprecise 
parametric  value  in  the  experiment.  Further  investigation  is  needed  for  this 
case. 

•  Figure  15:  This  corresponds  to  runs  4917  and  4918  for  a  40  m,  1.25  A,  1,000 

ll 


Hz  electric  source.  Data  from  the  quasi-static,  Sommerfeld  (i.e.  WSU),  FDTD 
and  HFSS  methods  for  the  ^’-component  of  the  electric  field  are  plotted  as  a 
function  of  radial  distance  when  the  EMA  is  located  at  a  depth  of  8.23  m.  The 
exact  same  observations  that  were  made  for  Figure  13  apply  to  this  case. 

•  Figure  16:  This  corresponds  to  run  1016  for  a  190  A-m',  10  Hz  magnetic 
source.  Data  from  the  Sommerfeld  (i.e.  WSU),  FDTD  and  HFSS  methods  for 
they-component  of  the  electric  field  are  plotted  as  a  function  of  radial  distance 
when  the  EMA  is  located  at  a  depth  of  8.23  m.  In  this  case  the  FDTD  and 
Sommerfeld  methods  predict  correctly  the  experimental  outcome.  HFSS 
predicts  the  correct  trend,  but  the  data  are  corrupted  by  random  spikes.  The 
cause  of  these  spikes  is  understood.  Since  a  magnetic  source  will  weakly 
excite  an  electric  field,  the  ability  to  compute  such  weak  fields  is  strongly 
dependent  on  the  HFSS  mesh.  Further  work  on  mesh  generation  issues  is 
needed. 

•  Figure  17:  This  also  corresponds  to  run  1016  fora  190  A-m2,  10  Hz  magnetic 
source.  Data  from  the  Sommerfeld  (i.e.  WSU),  FDTD  and  HFSS  methods  for 
the  z-component  of  the  magnetic  field  are  plotted  as  a  function  of  radial 
distance  when  the  EMA  is  located  at  a  depth  of  8.23  m.  In  this  case  the  FDTD 
and  Sommerfeld  methods  predict  the  same  outcome.  HFSS  also  predicts  the 
correct  outcome,  but  there  are  some  noticeable  deviations  from  the 
experimental  data.  All  three  methods  are  in  slight  disagreement  with  the 
experimental  data. 

•  Figure  18:  This  corresponds  to  run  1117  for  a  0.5  A-m,  1 00  Hz  electric  source 
(i.e.  both  electrodes  are  on  the  same  boat).  Data  from  the  quasi-static, 
Sommerfeld  (i.e.  WSU),  FDTD  and  HFSS  methods  for  the  .v-  and 
z-components  of  the  electric  field  are  plotted  as  a  function  of  radial  distance 
when  the  EMA  is  located  at  a  depth  of  295  m.  (Note:  Two  different  FDTD 
codes  were  written  for  this  scenario  to  gauge  their  relative  robustness;  the 
results  are  virtually  the  same.)  Here  we  see  that  the  data  from  all  five  methods 
generate  virtually  the  same  result.  However,  none  of  the  methods  correctly 
predict  the  dip  in  the  experimental  field  data.  The  cause  of  this  is  due  to 
positional  ambiguity  of  the  moving  source.  For  close  range  observations  (in 
this  case  10  meters  or  so),  the  source  could  have  moved  several  meters  during 
the  one  second  sampling  interval,  which  suggests  that  its  effective  static 
position  is  unknown.  Since  several  meters  is  a  large  fraction  of  10  meters,  this 
positional  ambiguity  is  large.  In  earlier  cases  the  observations  were  on  the 
order  of  hundreds  of  meters  and  a  positional  ambiguity  of  a  few  meters  is 
inconsequential.  We  also  note  the  quasi-static  method  compares  favorably  to 
the  experimental  data,  which  is  to  be  expected  for  such  close  range 
observations. 
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•  Figure  19:  This  corresponds  to  runs  5035  and  5038  for  a  50  A-m,  10  Hz 
electric  source  (i.e.  both  electrodes  are  on  the  same  boat).  Data  from  the 
quasi-static,  Sommerfeld  (i.e.  WSU),  FDTD  and  HFSS  methods  for  the  x-  and 
z-components  of  the  electric  field  are  plotted  as  a  function  of  radial  distance 
when  the  EMA  is  located  at  a  depth  of  295  m.  Since  the  dip  in  the  field  data 
occurs  around  200  meters,  the  positional  ambiguity  effect,  as  noted  previously, 
is  not  present.  Also,  we  see  that  the  FDTD  data  and  the  experimental  data 
agree  the  best,  although  the  other  methods  give  reasonable  results. 

Recall  that  the  experimental  data  is  gathered  in  sets  of  one  second  sample  bins.  Each 
bin  corresponds  to  an  approximate  location  of  the  boat.  The  actual  position  of  the  boat 
is  ambiguous  due  to  a  non-zero  boat  speed.  Consider  Figure  20,  which  shows  the 
effect  of  using  bin  sample  sizes  in  excess  of  one  second.  As  the  sampling  interval 
increases,  the  expected  outcome  of  smoother  data  is  also  seen.  For  long-range 
observations,  the  increase  in  sample  bin  sizes  does  not  have  an  adverse  effect  on  the 
data.  This  will  not  be  case  for  close-range  distances  since  a  10  second  bin  size  will 
correspond  to  a  20  to  30  meter  positional  ambiguity  of  the  source. 
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Scenario  #1  -  Skiff  towed  @  40m  - 1000  Hz  -  50  A*m 
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Figure  13 


Scenario  #2  -  Skiff  towed  @  40m  - 1 00  Hz  -  48  A*m 

Runs  3046,  3044,  3042 


2e-5 

1  e-5 

6e-6 
4e-6 

2e-6 

1  e-6 

6e-7 

4e1oo 


400 


500  600  700  800  900  1000 

Distance  (m) 


Figure  14 
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E-Field  Amplitude  (V/m) 


Scenario  #3  -  Skiff  -  1 000  Hz  -  50  A*m 
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Figure  15 


Scenario  #4  -  Magnetic  source  - 1 0  Hz  - 1 90  A*m; 
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Figure  16 


15 


E-Field  Amplitude  (V/m)  H-Field  Amplitude  (A/m) 


Scenario  #4  -  Magnetic  source  - 1 0  Hz  - 1 90  A*m‘ 
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Figure  17 


Scenario  #5  -  Electric  Hull  - 1 00  Hz  -  0.5  A*m 


Figure  18 
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E-Field  Amplitude  (V/m) 


Scenario  #6  -  Skiff  @  40m  - 1 0  Hz  -  50  A*m  -  EMA  @  295  m 
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Figure  19 
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Figure  20 
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Future  Work 


Based  on  the  previous  discussions  and  the  objectives  of  the  project,  the  following 
future  work  is  envisioned: 

•  HFSS  predicts  the  correct  result  once  a  frequency-dependent  calibration  factor 
is  specified.  We  surmise  that  this  factor  is  related  to  the  wave-number,  since 
the  factor  is  proportional  to  f1  2.  Work  is  needed  to  understand  why  this  factor 
is  needed  and  its  precise  relationship  to  frequency. 

•  There  are  two  ways  to  implement  a  source  in  HFSS;  each  way  has  its 
advantages.  We  need  to  investigate  whether  Hertzian  or  scattered  field  sources 
are  best. 

•  HFSS  requires  us  to  import  Lake  topographical  and  bathymetry  data.  Since  the 
Lake  floor,  in  particular,  is  smooth  in  some  places  and  spatially  varying  in 
other  places,  we  need  to  optimize  the  HFSS  grid  to  these  places  to  assure  fast 
simulations  and  accurate  data.  Mesh  optimization  will  be  done  with  a 
specialized  mesh  editor  from  Sandia  laboratories. 

•  HFSS  introduces  errors,  particularly  in  the  magnetic  field  calculation  at  ELF 
frequencies.  This  calculation  can  be  improved  by  mesh  refinement.  More 
research  is  needed  to  figure  out  how  the  mesh  should  be  refined. 

•  Boat  motion  errors  associated  with  source  position  ambiguity  are  made 
manifest,  particularly  close  to  EMA.  Further  work  is  needed  to  account  for  this 
effect  in  the  modeling  software. 

•  FDTD  results  tend  to  be  slow  on  single  processor  or  single  core  machines. 
Multi-core  FDTD  optimization  and  parallelization  of  the  FDTD  codes  will  be 
implemented. 

•  Mesh  truncation  schemes  for  FDTD  have  been  explored  in  Phase  1.  Further 
work  in  the  area  of  perfectly  matched  layers  is  still  needed. 

•  Ansys/Ansoft  has  a  product  called  “Maxwell,”  which  is  a  static  field  solver. 
Research  is  needed  to  see  if  this  code  can  predict  the  quasi-static  results. 

•  The  Sommerfeld  methodology  has  been  shown  to  predict  accurate  results.  One 
advantage  of  a  closed-form  solution  is  the  ability  to  predict  the  up-over-down 
effect.  This  effect  suggests  that  the  up-over-down  wave  is  less  attenuated  (with 
respect  to  the  direct  wave)  once  the  wave  escapes  into  the  air  medium.  This 
analysis  will  be  accomplished  in  Phase  Two. 

•  An  experiment  in  the  Fall  of  2009  and  the  Spring  of  2010  is  scheduled  to 
produce  more  validation  data. 

•  During  Phase  Three  user  manuals  and  documentation  will  be  written  for  each 
of  the  delivered  models  and  codes. 
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Appendix  A:  Sommerfeld  Half  Space  Method 


Robert  G  Olsen,  Zhi  Li 

School  of  Electrical  Engineering  &  Computer  Science,  Washington  State  University 

Formulation 

The  purpose  of  this  project  is  to  find  analytical  solutions  for  the  electric  and  magnetic 
fields  in  any  layer  of  a  3-layer  non-magnetic  media  model  when  a  vertical  or 
horizontal  magnetic  dipole  (VMD  or  HMD)  is  radiating  in  the  middle  layer.  The  three 
layers,  from  top  to  bottom,  are  assumed  to  be  air,  lake  water  and  lake  bottom.  And 
they  are  denoted  as  layer  #0,  #1  and  #2,  respectively.  The  illustrations  of  the  model 
are  shown  in  Fig.  1. 


,z 

,Z 

Air 

°b»  Po  Air 

to,  °0>  A) 

z  =  0 

z  =  o  y 

Water 

observation  Water 

observation 

VMD  it  z  -  -h 


point 


HMD  =*  z  =  -h 


point 


£i>  °i.  Md 

£i.o-i.Po 

777777777 

2  =  ^  V ////////// 

777777777 

Z  V.  d  '7777/777777 

Bottom 

h'  ^ 0 

Bottom 

£2.  ^2-  ft 0 

Vertical  magnetic  dipole  Horizontal  magnetic  dipole 

Fig.  1  Model  illustrations 

As  marked  in  the  figure,  £,  (/  =  0,  1  and  2)  is  the  permittivity,  //,  is  the  permeability 
and  a,  is  the  conductivity  of  the  ih  layer.  £l-£n£0,  where  sn  is  the  relative 

permittivity  and  so  is  the  permittivity  of  free  space.  For  non-magnetic  media,  /*/- 

Ho. 

An  integral  form  of  the  electric  vector  potential  is  used  to  derive  the  solutions  to  the 
fields.  For  VMD  and  HMD  cases  the  vector  potentials  can  be  respectively  written  as 
(assuming  the  time  variation  is  e?wt) 
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VMD 
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In  these  expressions  of  the  vector  potentials,  the  constants  kt  and  7,  and  the  variable  u, 
are  defined  as  (i  =  0,  1  and  2): 

k  jope’JdA 
1  An 
7, 2  =  ~or  ps, 


ui  -  yj(A2  +7(2) 


•  O', 


where  £,'  =  £,  —j—  is  the  complex  permittivity  of  the  ib  layer,  Re(y,)<0  and 

co 

ehR 

Re(j/J>0.  The  term  k] -  in  F]  for  VMD  and  F]  for  HMD  represents  the 

R 

radiation  field  by  the  dipole  in  an  infinite  homogeneous  medium  (infinite  lake  water 
for  this  case).  It  can  be  written  in  integral  form  as 


R 


k,  £  u^e'“,(:*h)AJ0 (Ap)dA  (z  +  h)  >  0 
Jfc,  £  u~' e“,u+h) AJ0(Ap)d A  (z  +  h)  <  0 


where  R  =  (xi+y2)v2. 


The  coefficient  functions//  through  /  and  gi  through  gv  can  be  determined  from  the 
tangential  components  of  the  electric  and  magnetic  fields  E  and  H ,  respectively,  at  the 
interfaces  (z  =  0  and  z  =  -d  planes),  They  are 
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where  DEN ,  =  (w0  —  m,)(m2  —  ;/,)e  —  (w0  +w,)(wl  +  u2)eu'J . 


With  these  coefficient  functions  being  determined,  the  analytical  solutions  to  the  E 
and  H  fields  can  be  derived  from  Maxwell's  equations.  The  complete  expressions  for 
E  and  H  fields  in  the  water  layer  are  listed  below. 
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VMD:  E  and  H  fields  in  water  (Medium  #1): 
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IdA 
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_[  X2u]'e~",(:+h)J0(Xp)dX  (z  +  h>  0) 
£  X3tt;'eu'(!+h)JQ(Xp)dX  ( z  +  h<  0) 


+  [[e^f2^f2]X%{Xp)dX 


HMD:  E  and  H  fields  in  water  ( Medium  #  1 ): 


flf 
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dFl_5F[ 

dy  dz 

£  -Xe"‘'{:+h)J0(Xp)dX  (. z  +  h>  0) 
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+  [  (e"“'rg2+e“lIg3) 

dF] 


+  J,  [-e~"'!f2  +  eu,:f2]"Mo('Lp)dX 
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y  el  dx 
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p 


X1  sin(0)cos(0V  X 
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e\  dx 
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u7leu'i:+h)  (z  +  h<  0) 


+  e~"':f2  +e“':f3  -ufi~K,:g2+ufivgi 


A  J0  (Ap)  sin2  (<fi)  h —  J,  (Ap)  cos(20) 
P 


A2dA 


( \ ,r'‘>-"'{z+h) 

♦rf 1  ■ 
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i/|-le",<*+*>  (z  +  /3<0) 
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CO/J£x 
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+e~^f2+e^f3 


y2  Ad 0{Ap)d  A 


dz 


dFv  +  dF[: 


dy  dz 


-r\K 


Id  At 
An 


l'-e'",(z+*)  (z  +  h>0) 


,«i(*+A) 


(z  +  /?<0) 


•  A2./,  (Ap)  sin(^)  JA“ 


-«,e-"':/2  +i/1e"'-/3  +  AV"zg2  +AV“-g3 


(/>  =  atan(x  I  y) 

With  the  analytical  solutions,  the  electric  and  magnetic  fields  in  the  water  layer  can  be 
easily  calculated.  Fig.  2  through  Fig.  5  show  some  examples  of  the  field  calculation. 
The  parameters  used  are  listed  in  the  table  below. 
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Air 

Water 

Bottom 

Relative  permittivity, 

ert(e,=e0  *„) 

1 

1 

1 

Conductivity,  a  (S/m) 

0 

0.018 

0.012 

Permeability,  n  (H/m) 

4jix  10’7 

4TIX10'7 

4rtxl077 

</(m) 

300 

ft  On) 

can  vary  from  0  to  300 

Dipole  moment  IdA  ( A  •  m  ) 

2500 

frequency /(Hz) 

1000 

Since  H-_  and  Hy  are  respectively  the  dominant  components  for  VMD  and  HMD  cases, 
the  results  for  only  these  two  components  are  plotted  here.  Fig.  2  and  Fig. 3  give  the 
magnitude  of  H:  for  the  VMD  case.  Given  the  measurement  equipment  available,  the 
minimum  detectable  magnetic  field  is  also  indicated  in  the  figure  (dotted  line).  From 
experimental  data.  |H|mjn  =  4x  10's  A/m  and  |E|mjn  =  lxl 0*^  V/m. 


Fig.  2  Magnitude  of  Hz  in  water  for  VMD  case:  VMD  is  5m  below  surface  (h  =  5m),  observation 

points  horizontally  placed  at  z  =  ~20m 
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Fig.  3  Magnitude  of  Hz  in  water  for  VMD  ease:  VMD  is  5m  below  surface  (h  =  5m),  observation 

points  vertically  placed  at  p  =  200m 

Fig.  4  and  Fig.  5  give  the  magnitude  of  Hy  for  the  HMD  ease: 
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Magnitude  of  H  (A/m) 


Fig.  4  Magnitude  of  Hy  in  water  for  HMD  ease:  HMD  is  5m  below  surfaee  (h  =  5m),  observation 
points  horizontally  placed  at  z  =  -20m,  O  =  rc/3 


Fig.  5  Magnitude  of  Hy  in  water  for  HMD  ease:  HMD  is  5m  below  surfaee  (h  -  5m),  observation 
points  vertically  placed  at  p  =  200m,  <t>  =  n! 3 
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Validation 


After  the  analytical  solutions  were  developed,  several  tests  were  conducted  to  validate 
the  results. 

VMD  Case: 


First,  the  formulas  for  the  field  were  checked  against  to  the  results  obtained  in  some 
references.  Reference  [1]  provides  formulas  for  the  calculation  of  magnetic  fields 
inside  a  semi-infinite  conducting  medium.  If  let  £2  -  £;and  02  =01,  our  3-layer  model 
reduces  to  a  2-layer  model,  which  is  analogous  to  the  model  used  in  [1].  Then  our 
equations  of  the  coefficient  functions  reduce  to  equation  (11)  and  (12)  in  [1],  A 
MATLAB  program  was  written  to  numerically  compute  our  field  formulas.  Again,  if 
the  model  is  reduced  to  a  2-layer  (air-water)  model,  our  simulation  result  matches 
with  Fig.  2  in  [1  ]. 


Fig.  6  Comparing  our  numerical  simulation  results  with  Fig.  2  in  [  1  ] 

To  facilitate  comparison,  a  factor  G  is  introduced  that  is  proportional  to  H:,  which  is 
also  used  in  [1]: 

//.  =  — ^-rG 
2n{z  +  h ) 

Likewise,  R 1  is  a  ratio  defined  as 
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where  z  is  the  depth  of  the  observation  point  and  h  is  the  depth  of  the  dipole  source.  If 
let  e o  =  c/and  an  =07,  the  3-layer  model  reduces  to  a  2-layer  (water-bottom)  model. 
The  same  simulation  results  for  the  factor  G  were  obtained  as  shown  in  Fig.  2  of  [I] 
when  the  dipole  source  was  put  at  the  water-bottom  interface. 

Second,  the  results  reported  here  were  also  compared  to  those  for  a 
homogeneous-medium  case.  For  instance,  if  the  VMD  source  is  placed  in  the  middle 
of  the  second  layer,  where  it  is  far  away  from  both  interfaces,  our  field  results  are  very 
close  to  that  of  a  quasi-static  case  in  which  a  dipole  radiates  in  the  infinite 
homogeneous  medium.  Figure  7  shows  an  example  for  this.  For  the  3-layer  case,  the 
dipole  is  150m  away  from  both  of  the  interfaces.  The  calculation  of  field  perfectly 
matches  with  that  for  the  quasi-static  case. 


h=-150;z  =  -150 


Fig.  7  VMD  simulation  result  compared  to  homogeneous  ease 

Third,  the  boundary  conditions  across  both  the  air-water  and  water-bottom  interfaces 
were  checked.  It  is  not  difficult  to  find  the  solutions  to  the  fields  in  the  first  and  third 
layer  of  the  medium  using  the  same  method  described  in  part  I.  Figures  8a  and  8b 
show  the  boundary  conditions  data  at  the  air-water  interface  for  the  VMD  case. 
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Magnitude  of  H  (A/m)  Magnitude  of  E  (V/m) 


h  =  -10  z  =  0 


(a)  E* 

h  =  -10  z  =  0 


(b)  Hp 

Fig.  8  Checking  tangential  E  &  H  fields  along  the  air-water  interface  for  VMD  case 
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The  continuity  of  tangential  fields  along  the  water-bottom  interface  was  also  validated. 
The  results  are  not  be  plotted  here. 

HMD  Case: 


For  HMD  case,  the  same  strategies  as  in  VMD  case  were  used  to  validate  the  results. 
First,  when  the  model  was  reduced  to  the  2-layer  case,  our  equations  for  the 
coefficient  functions  reduced  to  equations  (21)  through  (23)  in  [2],  Then  the 
simulations  shown  in  Fig.  2  in  [2]  were  repeated.  Our  numerical  simulations  (for 
fields  in  line  with  dipole,  refer  to  Fig.  2  in  [2])  of  the  factor  A  are  plotted  in  Fig.  9. 


<t>  =  0.5rc 


Fig.  9  Magnitude  of  factor  A  (proportional  to  Hy)  vs  relative  horizontal  distance  (0^0. Stc) 
The  factor  A  is  proportional  to  the  Hv  field.  In  [2],  it  is  defined  as 

A(yh)  =  ^-H^ 

V/  '  A  > 
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It  should  be  mentioned  that  our  results  are  not  exactly  the  same  as  those 
given  in  [2],  The  starting  points  and  the  final  values  match.  But  there  is  an 
overshoot  in  Fig.  9,  which  cannot  be  seen  in  Fig.  2  of  [2],  The  reason  for 
this  appears  to  be  that  in  [2]  a  first  order  approximation  of  integral  was 
used  to  calculate  the  field,  while  our  simulation  used  the  analytical 
solution,  which  can  give  more  accurate  calculation.  To  verify  that  ours  is 
correct,  the  boundary  conditions  for  this  specific  case  were  also  checked. 
From  our  numerical  data,  the  factor  A  (which  is  proportional  to  tangential 
H)  is  perfectly  continuous  across  the  interface,  as  shown  in  Fig.  10.  This 
provides  supporting  evidence  to  that  our  data  are  more  accurate  than  that 
in  [2].  The  tangential  E  fields  across  the  interface  were  also  checked  and 
the  similar  results  are  shown  in  Fig.  10. 


ItpI 


(a)  |//;|=0.l 
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(b)  |}'/, |=0.2 

Fig.  10:  Boundary  conditions  associated  with  factor  A  (0=0. 5 Jt) 


The  results  of  HMD  were  also  compared  to  those  for  the  quasi-static  case.  They  also 
match  with  each  other  very  well. 

Finally,  the  boundary  conditions  were  checked  for  both  tangential  electric  and 
magnetic  fields  across  air-water  and  water-bottom  interfaces.  The  results  for  the 
electric  field  at  the  water-bottom  interface  are  shown  in  Fig.  1 1 . 
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h  =  290:4)  =  0.25* 


5 

p(m) 


14 


x  10 


12 


(a)  Ex 


h  =  290;<J>  =  0.25* 


(b)  Ey 

Fig.  1 1 :  Tangential  E  Helds  along  the  water-bottom  interface  for  HMD  case 
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Appendix  B:  FDTD  Method 


Dennis  Sullivan 

Department  of  Electrical  and  Computer  Engineering,  University  of  Idaho 
Abstract 

The  major  effort  of  the  FDTD  group  for  the  past  year  has  been  to  extend  the  range  of 
FDTD  simulations  underwater  without  using  enormous  computer  resources.  It  has 
included  the  development  of  a  near-field,  far-field  formulation  to  allow  the  use  of 
large  cells  in  simulating  the  far  field,  and  the  development  of  an  interpolation  scheme 
to  improve  accuracy  when  large  cells  are  used  in  the  far  field.  It  has  also  included 
the  development  of  a  new  perfectly  matched  layer  to  reduce  the  size  of  the 
surrounding  computational  domain  A  reformulation  of  the  basic  FDTD  method  has 
been  introduced  to  improve  stability  during  extremely  long  computation  times. 

Introduction 

In  the  past  year,  the  FDTD  simulation  group  has  accomplished  the  following: 

1 .  The  reformulation  of  the  FDTD  code  for  greater  stability  at  low 
frequencies. 

2.  The  development  of  a  near-field/far-field  formulation  that  allows 
the  simulation  of  energy  sources  in  the  near  field  w  ith  relatively 
small  cells,  while  the  far  field  is  simulated  with  larger  cells  to 
extend  the  range  of  the  simulation. 

3.  The  development  of  a  perfectly  matched  layer  (PML)  for 
simulations  in  lossy  media  at  ELF  frequencies. 

4.  The  preliminary  investigation  of  a  high-resolution  interpolation 
method  to  improve  resolution  in  the  far  field. 

Accomplishments: 

Details  of  the  above  accomplishments  are  given  in  the  following  sections. 

1 .  Reformulation  of  the  FDTD  for  greater  stability  at  low  frequencies. 

The  formulation  of  the  Maxwell’s  equations  that  is  normally  used  in  FDTD 
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simulations  is  [1] 


(1-la) 


8D 

dt 


=  Vx// 


D(co)  =  s'  (co)  E  (co) 

dH 


Mo¬ 


di 


=  -V  x  E 


(1.1  b) 


(1.1c) 


This  formulation  allow  a  great  deal  of  flexibility  in  simulating  complicated  materials 

because  the  various  materials  are  specified  by  the  complex  dielectric  constants  (ft>) 

and  having  a  separate  equation  to  calculate  the  electric  field  E  from  the  flux  density  D 
allows  one  to  bring  to  bear  signal  processing  techniques  for  the  formulation  in  FDTD. 
The  simulation  of  a  lossy  material  leads  to  equations  of  the  form 


El 


D"  -  Sn 
At  ■  CT 


e.  +  - 


sn  =  s" 1  +  e:  . 


(1.2  a) 


(1.2  b) 


This  formulation  is  very  stable  at  radio  frequencies  and  above,  but  leads  to  potential 
instabilities  at  ELF  frequencies.  The  reason  is  Eq.  (1.2  b),  which  is  the  time  domain 
implementation  of  an  integration.  An  ELF  FDTD  simulation  is  over  a  fraction  of  a 
wavelength,  which  means  that  an  integration  like  Eq.  (1.2b)  is  constantly  increasing. 
Over  long  simulations  on  the  order  of  20,000  iterations,  this  has  occasionally  led  to 
instability. 

Alternatively,  we  begin  with  the  following  time-domain  Maxwell's  equations 

8E 


,_0  =  V  x  //  -  ctE 

r  0  dt 

dH  ^ 

un - =  -V  x  E . 

0  dt 


(1.3  a) 


(1.3  b) 


We  assume  the  cell  size  is  Arand  the  time  step  is  At.  The  £"*'  can  now  be 
calculated  from  [2,  3] 

E (,+l  (k)  -  ca(k) E *  (k)  +  cb{k)\_Hn;'n (k  +  1/2)-  Hn;V2  (k  - 1/2)],  (1 .4  a) 

where 

,  cb  =  ca- — — — r*  (l.^b) 

(ATAy) 

There  is  a  crucial  choice  that  was  made  here.  Usually,  the  Ex  term  next  to  the 


ca  = 


1  + 


At  a 


£ 
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conductivity  is  averaged  across  the  two  time  steps 


r  E-\k)~E:(k) 

(C  (*)+£(*)] 

i  ir;1 2 (k +  i/2)- Hnvtin(k-\/2)^ 

c?  ~r 

A/ 

j 

2  J 

<i 

uf 

which  would  lead  to  the  following  expression  for  ca : 


ca  = 


1- 


A/  •  o 

J 


1  + 


A  t-<7^ 

j 


(1.5) 


At  ELF  frequencies  in  lossy  media,  the  ca  of  Eq.  (1.5)  would  be  negative,  leading  to  a 
potentially  unstable  condition.  (The  implementation  of  Eq.  (1.3  b)  into  FDTD  is 
straight-forward  and  will  not  be  presented  here.) 


There  is  another  choice  that  leads  to  substantially  larger  time  steps,  and  therefore, 
substantially  faster  solutions.  Once  the  cell  size  Ax  is  chosen,  the  time  step  must  be 
chosen  to  satisfy  the  Courant  condition,  which  in  three  dimensions  is 

(1.6) 

where  cmix  is  usually  the  speed  of  light  in  a  vacuum.  The  complex  dielectric 
constant  is  calculated  by 

.  a 

£,  =£r  +~. - 

The  materials  of  interest  for  this  project  and  their  dielectric  properties  are  listed  in 
Table  one.  At  ELF  frequencies,  the  imaginary  part  of  the  dielectric  constant  will 
dominate  the  magnitude  for  all  the  materials  except  air.  Therefore,  increasing  the 
dielectric  constants  of  mud  or  metal  to  80  would  make  very  little  difference.  If  we 
assume  every  material  in  Table  one  has  a  real  dielectric  constant  of  80,  then 

cmx  =c0/\lso  and  the  time  step  is  almost  an  order  of  magnitude  greater.  (c0  is  the 

speed  of  light  in  a  vacuum.)  Even  though  air  is  one  of  the  materials  used  in  the 
simulations,  air  is  a  boundary  medium  in  this  project.  It  presents  almost  perfect 
reflection  to  an  electromagnetic  signal  in  water,  even  if  the  higher  dielectric  constant 
is  used. 


A  t< 


Ax 


vr; 


Table  1.  The  properties  of  the  materials  used  in  the 
simulations  described  in  this  paper  [4], 


Material  er  cr  (S / m) 

Air  1  0 

Lake  water  80  0.0 1 8 

Mud  40  0.002 

Metal  1  107 
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2.  The  development  of  a  near-field/far-field  formulation 

By  utilizing  the  equivalence  principle  [5],  a  method  has  been  developed  to  model  the 
radiating  sources  with  a  relatively  high  resolution,  and  to  model  the  far  field  with 
larger  cells  to  extend  the  potential  problem  space.  Two,  three-dimensional  FDTD 
problem  spaces  are  utilized  to  implement  the  near-to-far  field  transformation  (Fig.  2.1) 
[6,  7].  A  problem  space  with  a  relatively  small  cells  size  (1  m‘)  is  used  to  model  the 
source  (Fig.  2.1a).  Another  problem  space  with  larger  cells  (13  m3)  is  used  to  model 
the  far  field  (Fig.  2.1b).  The  ratio  of  13  to  1  between  far  and  near  field  cells  size  is 
somewhat  arbitrary.  For  the  simulations  in  this  paper,  it  resulted  in  near  and  far  field 
problem  spaces  of  about  the  same  size.  Each  problem  space  is  surrounded  by  a 
perfectly  matched  layer.  (A  new  PML  for  ELF  frequencies  and  lossy  media  has  been 
developed  and  is  discussed  in  section  3.)  Each  problem  space  contains  a 
three-dimensional  transfer  surface  where  the  equivalence  principle  is  implemented. 
The  tangential  fields  calculated  on  the  transfer  surface  in  the  near  field  are  impressed 
on  the  transfer  surface  in  the  far  field  to  form  the  far  field  source.  This  is  illustrated 
in  Fig.  2.2. 


Transfer 

surface 


Figure  2.1.  Two  problem  spaces  are  used  in  the  FDTD  simulation. 


-6 

x  10 

1, 


(a)  (b) 

Figure  2.2.  (a)  Near-field  mesh  plot,  (b)  Far-field  mesh  plot. 
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Validation 


In  this  section,  the  results  of  the  near  to  far  field  transformation  are  compared  to 
analytic  results  calculated  using  Sommerfeld’s  half-space  (SHS)  problem. 
Sommerfeld’s  half  space  problem  calculates  the  resulting  fields  from  an  oscillating 
dipole  near  a  plane  interface  separating  two  homogeneous  half-space  regions,  as 
illustrated  in  Fig.  2.3.  This  method  was  well  described  in  Appendix  A  and  in  the 
literature  [8]  and  will  not  be  repeated  here.  In  Fig.  2.3,  the  upper  layer  is  air,  the 
middle  layer  is  water  and  the  lower  level  is  mud.  The  water  layer  in  the  middle  is 
300  meters  thick.  The  monitor  lines  represent  the  plaees  where  comparisons 
between  the  methods  will  be  made. 


Figure  2.3.  Diagram  of  the  comparison  between  the  FDTD  near-to-far  field 
transformation  and  the  SHS  method.  The  souree  is  1  meter  below  the 
air-water  interface. 

Comparisons  at  200  meters  and  1000  meters  are  shown  in  Fig.  2.4a  and  Fig.  2.4b, 
respectively.  The  amplitudes  in  each  figure  are  calculated  by  the  method  of 
two-equations,  two-unknowns  [3].  The  horizontal  coordinate  is  the  distance  to  the 
air/water  surface  and  the  vertical  coordinate  is  the  magnitude  of  the  field.  The 
discrete  symbols  represent  the  FDTD  calculations  and  the  different  kinds  of  lines 
represent  the  calculations  by  SHS  method.  Clearly,  the  results  of  the  comparisons 
are  very  good  in  all  cases. 


Distance =200  m 


Distance=1000  m 


(a) 


FDTD  10HZ 
+  FDTD  100Hz 
O  FDTD  1kHz 
SHS  10HZ 
-  -  SHS  100Hz 
SHS  1kHz 


t.~»r  rr  -r  ** ± ■ 
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Figure  2.4.  Comparisons  of  the  FDTD  simulations  and  the  SHS  calculations  for  the  cell 
size  ratio  of  10.  The  souree  is  near  the  upper  surface  of  the  water  layer,  (a)  200  meter 
ease,  (b)  1000  meter  case. 
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This  section  illustrates  the  use  of  FDTD  with  a  near  field,  far  field  transformation  in 
simulating  a  more  realistic  case.  One  of  the  goals  of  the  simulation  is  to  verify  the 
accuracy  of  the  modeling  with  measured  data  that  will  be  made  in  a  lake.  Because 
the  lake  bed  is  not  flat,  analytical  approaches  cannot  be  used  to  solve  the  problem. 

The  antennas  that  will  be  used  in  the  transmission  are  four  meter  by  four  meter 
rectangular  current  loops.  Figure  2.5a  illustrates  how  this  is  simulated  in  the  XY 
plane  in  the  FDTD  space.  Metal  can  be  simulated  by  insuring  that  an  E  field  is  zero 
at  a  particular  point  in  the  space.  Therefore,  using  cells  that  are  one  meter  cubed,  the 

metal  loop  antenna  is  simulated  at  the  corresponding  Ex  or  Ev  positions,  as  shown 

in  Fig.  5a.  Since  the  radius  of  the  w  ire  of  the  antennas  is  considerably  less  than  the  1 
meter  cells  size,  the  thin  rod  approximation  [9]  is  used  to  model  the  wire  at  these 
positions. 

In  FDTD,  a  current  cannot  be  simulated  directly,  but  it  can  be  simulated  indirectly  by 
using  Ampere's  circuit  Law  [5]  and  specify  ing  the  surrounding  H  fields: 

/  =  cf—  H  •  c/I .  (2.1) 

JEo 


By  impressing  a  hard  source  on  one  of  the  Ev  fields,  a  value  is  induced  on  the 
surrounding  H  fields,  as  shown  in  Fig.  2.5b.  This  results  in  a  current,  via  Eq.  (2.1). 
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Figure  2.5.  The  simulation  of  a  four  meter  by  four  meter  current  loop  in  FDTD. 

A  model  of  the  lake  bed  is  created  for  the  far  field  domain  (Fig.  2.6).  The  cells  are 
ten  meters  cubed.  The  shape  of  the  lower  surface  of  the  water  layer  shows  a 
complex  geometry  structure  similar  to  a  real  lake  bed. 

Figure  2.7  shows  the  results  of  the  near  to  far  field  simulation  using  the  current  loop 
source  in  the  near  field  and  the  lake  bed  in  the  far  field.  Results  are  given  for  two 
different  places,  200  meters  and  1000  meters  from  the  source,  and  at  three  different 
frequencies,  10,  100.  and  1000  Hz.  These  simulations  were  done  on  an  HP  DL140 
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GE  Quad  Core  and  required  about  6  hours.  Both  the  near  and  far  field  problem 
spaces  were  120  cells  cubed. 


Figure  2.6.  The  lake  bed  that  is  simulated  in  the  far  field.  The  cells  are  10  meters 
cubed.  The  souree  is  the  eurrent  loop  located  one  meter  below  the  surface. 


Distai>ce=200  m 


Distanced 000  m 


(a)  The  results  at  200  meters 


(b)  The  results  at  1000  meters. 


Figure  2.7.  Results  of  the  simulation  illustrated  in  Fig.  2.6. 

3.  The  development  of  a  perfectly  matched  layer  (PML)  for  lossy  media  at  ELF 

frequencies 

Berenger  [10]  assumed  that  any  plane  wave  propagating  in  the  direction  d  near  the 
PML  could  be  broken  up  into  the  part  traveling  perpendicular  to  the  PME,  c/±,  and 

the  part  traveling  parallel,  d  (Fig.  3.1).  The  two  conditions  for  the  PML  are 

1.  It  must  have  the  same  impedance  as  free  space  and  not  present  a  loss  to  the 
wave  traveling  parallel  to  the  interface  [Eq.  (3.1)]. 
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2.  It  must  increase  the  artificial  electric  and  magnetic  conductivities  such  that  the 
impedance  still  matches  that  of  the  free  space. 

Both  of  these  conditions  are  met  by  increasing  the  electric  and  magnetic 
conductivities  in  the  PML  such  that 


Note  that  this  impedance  is  a  real  number. 

Berenger  implemented  Eq.  (3.1)  into  the  FDTD  formulation  by  a  split-step 
formulation  that  broke  each  electric  and  magnetic  field  up  into  two  components. 
Most  applications  assume  that  the  background  medium  in  the  main  problem  space  is 
free  space. 


PML  Medium 


Free  space 

V**0 


Figure  3.1.  The  PML  is  implemented  by  assuming  any  propagating  wave  can  be  broken  up 
into  a  part  that  is  perpendicular  to  the  PML  interface  and  a  part  that  is  parallel  to  it. 

The  PML  in  a  lossy  medium  at  ELF  frequencies 

When  the  background  medium  is  lake  water  and  the  frequencies  are  in  the  ELF  region, 
the  situation  is  different  [10,  11].  Lake  water  has  a  dielectric  constant  of  80  and  a 
conductivity  of  0.018  S/m.  For  lake  water  at  1  kHz,  the  complex  dielectric  constant 
is 


0.018 


=  80  -  y'3.24  x  1 0s  =  -y‘3.24  x  1 05. 
Therefore,  the  impedance  is 


(3.2) 
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Vtr=\nH\Z45°. 


Notice  that  because  the  loss  term  dominates,  the  impedance  is  45  degrees  in  phase. 
The  impedance  of  the  PML  material  must  remain  at  this  value,  but  at  the  same  time 
increase  the  loss  further  as  it  goes  perpendicular  into  the  PML.  This  can  be 
accomplished  by  adding  a  factor  s  to  the  conductivity  and  the  permeability 


n„, 


serw 
\  j(oe o 


(3.3) 


This  causes  the  PML  medium  to  absorb  outgoing  waves  faster  than  the  water  medium, 
but  also  avoids  reflections  from  the  PML  medium. 


Implementation  into  FDTD 

We  begin  by  looking  at  the  implementation  of/s,  and  II  propagating  in  the  y 
direction  perpendicular  to  the  PML: 


dE„  dH 

e...Sn— — +<r  £  = — L 


dt 


dy 


dH,  dE 


Note  that  sH  =  80  and  o  =  0.018  for  water.  The  first  step  to  convert  the  equation 

to  the  FDTD  formulation  is  to  approximate  all  partial  derivatives  as  difference 
equations.  The  spatial  derivatives  are  not  as  important  in  this  discussion,  so  we  write 

8H:  ~  del  _  H 
dy  Ax 

Taking  the  usual  finite-difference  approximation  to  the  time  derivatives  gives 

\  del _H 

t,  —  9 

V 


and  solving  for  the  updated  version  of  E"*' 

A t  del_H 
e„£0  Ax 

From  this  we  can  develop  the  FDTD  equations: 


E"*'  -  E"  + 


A  t-a„ 


£  £ » 


=  ca  •  E"y  +  cb  •  del _  H , 


(3.4  a) 
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ca  = 


1 


(3.4  b) 
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At  •  cr 
+  — 

V 


del  H' 


A  similar  procedure  for  the  //  field  gives 

H"*'n  =  H"~U2  +  db  ■  del  _  E , 

where 


(3.5  a) 


db  =  -^—.  (3.5  b) 

//0-Ax 

Recall  that  to  implement  the  PML  we  increase  the  conductivity  and  the  permeability 
at  the  same  rate  by  a  constant  parameter  s  as  shown  in  Eq.  (3.3).  This  is  easy  to  do 
for  the  permeability.  We  add  another  parameter  jy  to  Eq.  (5  a) 

HT'n  =  H"  1/2  +  fy{j)-db-del_E.  (3.6) 

Instead  of  changing  db  to  A//(2-/v0-Ax) ,  for  instance,  we  simply  set  fy(j)  =  0.5  . 


However,  increasing  <jw  is  not  so  straight-forward.  Notice  that  it  involves  both  ca 

and  cb.  Rather  that  recalculate  the  parameters  for  each  increase,  we  define  another 
parameter 


.  At  ■  a,  N 

I  + 


gy{j)  = 


o  7 


|  +  A,.s.^ 


(3.7) 


B, 


o  / 


The  factor  is  added  to  Hq.  (3.4  a): 

Enx;'  =  gy  (J)  ■  ca  ■  El  +  gy(j)  • cb-del  _H . 
Note  that  it  is  necessary  to  split  the  E  fields,  but  not  the  H  fields. 


(3.8) 


Results 

In  this  section  we  illustrate  the  effectiveness  of  the  lossy  medium  PML.  We  will 
start  with  the  problem  space  illustrated  in  Fig.  3.2a,  which  is  80  cells  cubed.  Each 
cell  is  25  meters  cubed.  The  source  is  a  single-cell  electric  dipole.  After  4000  time 
steps,  the  amplitude  is  calculated  along  a  monitor  line  as  shown  in  Fig.  (3.2  b). 
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Figure  3.2.  (a)  A  dipole  source  is  located  in  the  801  problem  space.  Once  steady  state 

has  been  reached,  the  amplitude  is  determined  at  a  transverse  line  15  cells  from  the  source. 

The  cells  are  25  m*.  (b)  The  amplitude  at  the  monitor  line  after  4000  time  steps. 

The  simulation  is  then  repeated  for  the  truncated  problem  space  shown  in  Fig.  3a 
where  a  four-cell  PML  has  been  added.  In  this  simulation  the  right  wall  has  been 
moved  in  to  within  ten  cells  of  the  source.  The  results  are  plotted  in  Fig.  3b  along 
with  the  results  of  the  previous  simulation.  For  comparison,  Fig.  3c  is  the  same 
simulation  with  no  PML  on  the  truncated  wall. 

In  one  final  simulation,  the  problem  space  is  reduced  to  60  x  20  x  20  cells,  as 
illustrated  in  Fig.  3.4a.  The  results  are  shown  in  Fig.  3.4b  where  the  results  of  the 
original  simulation  of  Fig.  2  are  presented  for  comparison.  The  amplitudes  on  the 
monitor  line  within  five  cells  of  the  center  are  identical.  If  these  values  represent  the 
information  of  interest,  then  the  addition  of  the  PML  means  a  problem  space  of  60  x 
20  x  20  =  24.000  cells  instead  of  the  original  803  =  510,000  cells.  Fig.  3.5c  repeats 
the  simulation  with  no  PML  on  the  boundaries. 

Because  this  PML  was  developed  for  lossy  media,  it  is  not  particularly  effective  in 
free  space  or  other  lossless  media.  Most  of  the  simulations  for  this  project  include 
the  water  surface  and  therefore,  free  space  is  part  of  the  problem.  What  remains  to 
be  developed  is  an  interface  mechanism  so  the  lossy  PML  can  be  used  underwater  and 
the  normal  PML  used  above  water. 
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Figure  3.3.  (a)  The  problem  spaee  is  truncated  to  ten  cells  to  the  left  of  the  source.  A 

four-cell  lossy  PML  has  been  added  to  each  boundary,  (b)  The  solid  line  is  the  amplitude  for 
the  simulation  in  (a),  while  the  dashed  line  is  from  the  80  eell  monitor  line  of  Fig.  3.2.  (c) 

The  same  simulation  with  no  PML. 
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Figure  3.4.  (a)  A  simulation  similar  to  Fig.  (3.2)  hut  with  the  problem  space  truncated  to  60  x  20  x 

20  cells,  (b)  The  solid  line  is  the  amplitude  for  the  smaller  problem  space  while  the  dashed  line  is  the 
larger  problem  space  of  Fig.  3.2.  (e)  The  same  comparison  when  a  60  x  20  x20  problem  spaee  with 
no  PML  is  used. 

4.  The  development  of  a  high  resolution  interpolation  method 

Section  2  described  the  near-field,  far-field  transformation  that  enables  us  to  use 
relatively  large  cells  to  extend  the  computational  domain.  However,  there  may  be 
times  when  relatively  high  resolution  is  needed  in  the  far  field.  Figure  4.1  illustrates 
such  a  possibility  where  an  indentation  in  the  sea  floor  would  be  lost  if  cells  on  the 
order  of  1 0  meters  cubed  are  used. 
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Figure  4.1.  Typical  contours  in  shallow  water.  Cell  sizes  on  the  order  often  meters  cubed 
are  used  in  the  far  field.  This  could  result  in  substantial  errors  when  modeling  shallow  water 
where  mines  are  often  planted. 

One  solution  is  a  correction  to  the  FDTD  fields  that  is  made  after  the  core  FDTD 
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simulation  is  finished  [13].  In  the  following  discussion  Zswill  represent  the  true 
electric  field  around  a  boundary,  and  E  represent  the  averaged  FDTD  electric  field 
value  (Fig.  4.2).  At  an  arbitrary  boundary,  the  electric  fields  can  be  represented  as 
the  sum  of  the  tangential  and  normal  components 

E  =  E,  +  En  (4.1) 

E  =  E,  +  En  (4.2) 


Figure  4.2.  FDTD  only  calculates  values  at  discrete  points,  such  as  Ez(i,j,k).  Values  at 
other  positions  can  be  approximated  by  interpolating 

Two  assumptions  are  made:  The  tangential  components  of  the  true  and  the  FDTD 
values  are  equal 


and  the  normal  flux  densities  are  equal 


E,  =  E,, 


D  =  D  . 

n  n 


Equation  (4.4)  leads  to  the  following  relationship  between  the  E  fields: 

eEn  =  , 


(4.3) 

(4.4) 

(4.5) 


where  e  is  the  true  dielectric  constant  at  that  point  and  e  is  the  averaged  dielectric 
constant  usually  used  in  the  FDTD  formulation.  From  Eq.  (4.3)  and  (4.5),  it  can  be 
shown  that  the  true  electric  field  can  be  calculated  by  the  averaged  FDTD  E  field  with 
a  correction  term  added.  Specifically,  in  the  z  direction  the  E  field  is 


E=E.+ 


£ 

—  -1 

ve 


(4.6) 


/ 


where  E,n  is  calculated  by 


E„  =  «.{('»  •£)'»]• 


(4.7) 


In  Eq.  (4.7)  ft  is  the  gradient  of  the  dielectric  constant  and  ft.  is  the  z  component  of 
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h .  The  correction  is  applied  after  the  FDTD  simulation  is  finished.  A  more 
detailed  explanation  is  in  [14]. 

Verification  of  the  Method 

In  order  to  verify  the  accuracy  of  the  interpolation  scheme  described  in  the  previous 
section,  we  used  an  analytic  method  based  on  Legendre  polynomials  [4].  This 
method  calculates  the  values  of  the  electric  field  inside  a  layered  sphere  in  an  electric 
field.  Table  one  is  a  list  of  the  materials  of  interest,  along  with  their  dielectric 
properties.  In  particular,  we  will  look  at  the  values  on  an  axis  at  forty-five  degrees  in 
the  YZ  plane  because  it  is  at  these  slanted  angles  that  the  largest  FDTD  error  occurs. 


Figure  4.3.  A  layered  sphere  is  used  to  provide  an  analytic  cheek  to  the  FDTD  data.  The 
background  medium  is  water;  the  sphere  is  half  mud,  half  water.  The  narrow  layer  in  the 
middle  is  either  mud  or  water. 

We  show  the  results  for  an  FDTD  simulation  using  10  m  cubed  cells  at  1  kHz  in  a 
water  medium  both  before  and  after  the  correction  is  applied.  The  sphere  in  Fig.  4.3 
has  the  characteristics  of  mud.  The  thin  layer  is  Vi  a  cell  wide  and  has  the 
characteristics  of  water.  Figure  4.4  displays  the  results.  FDTD  alone  sees  the  basic 
features  but  tends  to  average  the  magnitudes  out.  After  the  interpolation,  a  much 
better  agreement  is  attained,  even  though  the  layer  is  less  that  one  cell. 

Since  the  interpolation  takes  place  only  after  the  FDTD  simulation,  it  adds  nothing  to 
the  computation  time  and  very  little  to  the  needed  computer  resources.  This 
approach  provides  more  accuracy  when  large  cells  on  the  order  of  ten  meters  are  used 
for  long  range,  underwater  FDTD  simulation. 
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Figure  4.4.  Results  for  a  1/4  cell  (2.5  m)  layer  of  water.  The  straight  line  is  the  analytic 
values,  the  x's  are  from  and  FDTD  simulation  before  the  correction,  and  the  o‘s  are  the 
values  after  the  correction. 

Conclusion 

The  effort  over  the  past  year  has  concentrated  on  maximizing  the  effectiveness  of  the 
FDTD  simulation  over  long  distances.  Towards  this  end,  a  more  stable  formulation 
of  the  FDTD  algorithm  was  introduced  to  insure  stability  during  very  long  simulations. 
A  near-field,  far-field  transformation  was  implemented  to  allow  radiating  sources  to 
be  simulating  with  relatively  small  cells  to  maintain  accuracy,  while  the  far  field  is 
simulated  with  larger  cells  to  increase  he  simulation  distance.  A  perfectly  matched 
layer  was  developed  to  work  more  effectively  for  highly  lossy  materials  at  ELF 
frequencies.  This  decreases  the  need  for  a  large  surrounding  computational  space. 
Finally,  a  correction  to  the  FDTD  method  has  been  implemented  that  allows  the 
determination  of  the  fields  at  a  higher  resolution  than  the  FDTD  cells. 

One  of  the  tasks  that  remain  is  the  merging  of  the  new  underwater  PML  with  the  old 
PML  for  free  space.  The  underwater  PML  is  not  effective  in  free  space.  While  the 
propagation  of  ELF  signals  in  free  space  is  not  a  major  interest  for  this  project,  it 
could  play  a  substantial  role  in  maintaining  accuracy  of  the  underwater  simulations. 
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Appendix  C:  HFSS  Method 

Chris  Johnson  and  Christopher  Wagner 
Department  of  Electrical  and  Computer  Engineering,  University  of  Idaho 


Introduction 

Ansoft’s  High  Frequency  Structural  Simulator  (HFSS)  is  a  commercial 
frequency-domain,  finite-element  electromagnetic  fields  solver.  The  purpose  of  this 
particular  investigation  is  to  ascertain  whether  or  not  HFSS  is  applicable  for 
simulating  ELF  electromagnetic  signals.  Towards  that  end,  discussions  are  provided 
below  that  address  how  HFSS  should  be  used  and  validation  data  are  also  provided  to 
substantiate  the  efficacy  of  HFSS  for  this  application.  It  is  assumed  in  this 
discussion  that  the  reader  is  already  familiar  with  the  use  and  operation  of  HFSS;  this 
report  only  focuses  on  those  features  of  HFSS  that  are  unique  to  ELF  signal 
generation  and  data  processing.  Version  1 1.1. 3.1  was  used  in  the  following  discussion. 

Using  HFSS 

Hertzian  Sources 

The  Hertzian  dipole  is  a  source  option  available  under  the  excitations  menu  in  HFSS. 
Per  equivalence  theory,  this  source  is  capable  of  producing  the  fields  of  an  electric  or 
magnetic  dipole.  To  implement  this  source,  the  following  needs  to  be  executed: 

1 .  Select  “Excitation/Assign/Incident  Wave/Hertzian-Dipole  Wave" 

2.  Specify  the  location  of  the  source  in  the  model. 

3.  Select  the  dipole  moment  and  dipole  direction  “Cartesian/Spherical  Vector 
Setup”  tab. 

4.  Select  the  “Hertzian-Dipole  Wave  Options”  tab. 

5.  The  next  tab  has  two  options:  the  radius  of  surrounding  sphere  and  type  of 
dipole. 

The  first  option  needs  a  bit  of  clarification.  It  allows  the  user  to  specify  a 
mathematical  spherical  surface  that  surrounds  the  dipole.  On  this  surface  HFSS 
impresses  a  set  of  equivalent  currents  that  produce  the  same  field  as  that  of  a  dipole. 
Although  equivalence  theory  does  not  restrict  the  size  of  this  sphere,  the  size  does 
impact  the  simulation  efficiency.  A  small  sphere  requires  smaller  tetrahedrons  and  for 
a  given  domain  size,  this  implies  more  tetrahedrons  and  longer  simulations.  Larger 
spheres  imply  just  the  opposite,  but  preclude  examination  of  the  fields  close  to  the 
dipole,  since  one  cannot  ascertain  field  information  within  the  sphere.  Thus  a  tradeoff 
between  spherical  size,  computational  efficiency  and  computational  accuracy  is 


54 


required. 


For  ELF  applications  we  found  that  a  sphere  of  10  meters  does  not  cause  a  substantial 
increase  in  computation  time  and  allowed  sufficiently  close  observations  for  most 
applications.  If  a  study  of  the  effects  within  the  10  meters  is  necessary  then  a  second 
model  can  be  created  with  a  much  reduced  domain  so  that  the  computation  time  does 
not  become  extremely  burdensome.  Note:  HFSS  places  restrictions  on  the  size  of  the 
objects  and  domains  with  respect  to  the  default  units  of  the  model.  For  example, 
HFSS  does  not  allow  the  use  of  millimeter  sized  objects  when  dealing  with  kilometer 
sized  domain  and  vice  versa.  The  default  model  units  can  be  changed  under 
“Modeler/Units”  in  the  main  menu.  It  should  be  set  to  “km”  for  most  ELF 
applications. 

Edit  Source  menu 

HFSS  has  an  option  that  allows  the  user  to  change  how  the  fields  are  calculated 
during  the  post  processing  step.  Under  the  tab  HFSS\Fields\Edit  Source  there  is  an 
option  for  total  field  or  incident  field  calculations.  This  has  no  effect  on  how  the  fields 
are  calculated  during  the  solution  process  but  does  affect  how  the  information  is 
presented.  For  example,  the  scattered  fields  formulation  is  recommended  when 
dealing  with  plane  waves,  but  it  can  lead  to  errors  when  the  total  field  is  weak,  which 
suggests  that  the  incident  and  scattered  field  are  roughly  similar.  That  is,  the 
subtraction  of  the  incident  field  from  the  total  field  yields  a  small  number  that  is 
susceptible  to  numerical  noise.  This  error  is  particularly  pronounced  in  a  vacuum 
where  the  total  is  zero.  The  results  in  this  paper  were  all  done  with  the  total  fields 
option  selected. 

PEC  Block  Excitation 

The  Hertzian  sources  in  version  11.1.3  are  not  implementable  when  the  source  is 
embedded  in  lossy  media.  For  this  reason  there  was  considerable  doubt  early  on  about 
whether  HFSS  was  a  viable  tool  for  this  project.  Upon  further  investigation  and 
consultation  with  HFSS  technical  support  a  solution  to  this  problem  was  found.  It  was 
discovered  that  an  equivalent  source  can  be  constructed  using  a  perfectly  conducting 
block  of  any  shape.  This  is  accomplished  by  implementing  the  following  steps: 

1.  Place  a  perfect  electric  conductor  (PEC)  volume  where  the  source  is  to  be 
placed. 

2.  Embed  a  Hertzian  electric  dipole  or  magnetic  dipole  in  the  block. 

3.  Place  a  radiation  boundary  that  encompasses  the  block  using  the  “Incident 
Field"  specification. 

4.  Make  sure  the  surrounding  sphere  of  the  Herztian  source  is  completely 
contained  within  the  PEC  block. 
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We  surmise  that  this  approach  is  in  effect  using  a  scattered  field  formulation  by 
impressing  an  incident  field  on  the  block  and  allowing  the  block  to  scatter  that  field. 
The  scattered  field  is  equivalent  to  that  of  a  dipole.  This  approach  does  cause  a  scaling 
factor  to  be  introduced  in  the  data  that  is  still  under  investigation. 

Boundary  Conditions 

As  with  any  full  domain  solver  that  is  simulating  fields  within  an  open  domain,  care 
must  be  taken  to  truncate  that  domain  with  a  radiation  boundary  condition  that 
emulates  an  infinite  space.  HFSS  allows  the  user  to  choose  between  a  low  accurate 
radiation  boundary  condition  or  a  highly  accurate  perfectly  matched  layer  (PML).  The 
radiation  boundary  and  the  PML  boundary  are  both  available  under  the  “Boundaries” 
menu.  We  have  found  that  both  conditions  work  quite  well  for  ELF  applications.  Due 
to  the  simplicity,  we  typically  use  a  radiation  boundary  and  have  found  that  default 
settings  are  sufficient. 

Global  Materials 

The  global  material  setting  (typically  a  vacuum)  is  the  default  material  when  no 
material  is  selected  by  the  user.  However,  the  global  material  setting  also  impacts  the 
quality  of  the  radiation  boundary,  particularly  when  observation  points  are  chosen  that 
are  close  to  the  radiation  boundary.  HFSS  recommends  that  the  global  material  setting 
be  that  which  is  specified  for  the  interior.  For  example,  in  the  case  of  an  antenna 
radiating  in  a  vacuum,  the  global  material  setting  should  be  a  vacuum;  if  it  is  radiating 
in  water,  then  the  global  material  setting  should  be  water,  etc.  However,  in  a 
multilayered  or  inhomogeneous  material,  there  is  no  one  global  material  setting  that  is 
optimal.  To  obtain  the  best  solution,  however,  HFSS  recommends  that  the  global 
material  setting  be  that  of  the  material  that  surrounds  the  observation  point.  As  an 
example,  consider  a  three-layered  region  of  air,  water  and  mud.  If  the  observation 
point  is  in  air  and  close  to  the  radiation  boundary,  a  global  material  setting  of  air 
should  be  used;  if  the  observation  point  is  in  mud  and  close  to  the  radiation  boundary, 
a  global  material  setting  of  mud  should  be  used,  and  so  forth.  To  set  the  global 
material: 

1.  Select  the  model 

2.  Select  “HFSS/Boundaries/Edit  Global  Material  Environment” 

Mesh  SurfacesWolumes 

A  high-quality  mesh  is  critical  in  obtaining  good  results  from  HFSS.  A  sophisticated 
meshing  tool  is  embedded  within  HFSS  that  automatically  meshes  the  problem  space 
based  on  the  input  from  the  solid  modeling  editor.  However,  as  with  any  meshing  tool, 
the  resulting  mesh  is  not  always  optimal.  For  this  reason  it  is  sometimes  critical  to 
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guide  the  meshing  tool  by  adding  non-material  surfaces  that  assist  the  tool  to  place 
more  (or  less)  tetrahedrons  within  a  given  area  in  order  to  improve  accuracy  (i.e.  more 
tetrahedrons)  or  to  increase  computational  speed  (i.e.  less  tetrahedrons).  For  example, 
to  improve  the  mesh  along  a  field  monitoring  line,  a  rectangular  sheet  can  be  placed 
such  that  the  edge  of  the  rectangle  is  on  the  line  of  measurement.  Then  by  performing 
a  mesh  operation  on  the  rectangle,  the  user  can  increase  the  number  of  tetrahedrons  in 
that  localized  area  without  significantly  increasing  them  throughout  the  entire  volume. 
These  mesh  operations  are  available  in  the  project  manager  menu. 

Exporting  Field  Data 

Exporting  field  data  is  accomplished  using  the  “Fields  Calculator”  in  HFSS.  It  is 
located  under  HFSS\Fields\Calculator  within  the  main  menu.  Although  numerous 
options  are  available  in  the  Fields  Calculator,  only  a  few  are  relevant  for  this 
application. 

In  the  upper  right  corner  of  the  Fields  Calculator  is  where  the  specific  solution  set  of 
interest  is  specified.  Take  care  to  make  sure  that  each  parameter  is  set  properly  to 
ensure  a  valid  simulation,  per  the  procedure  below: 

1.  Generate  a  solution  using  the  “Analysis”  option  under  the  project  manager 
tree. 

2.  The  “Solution”  drop  down  menu  has  all  the  different  solutions  that  were 
generated  by  HFSS 

3.  After  the  solution  is  selected  then  any  field  type,  frequency  data,  or  phase  data 
that  has  a  saved  solution  can  be  selected. 

4.  If  a  parametric  sweep  was  done  then  the  different  parameters  can  be  selected 
with  the  “Change  Variable  Value”  button. 

After  specifying  the  solution  of  interest,  the  values  of  interest  need  to  be  specified. 

5.  Select  the  “Quantity”  button.  When  it  is  selected  a  drop  down  menu  will 
appear  and  E  and  //  values  can  be  selected.  When  either  is  selected  an  entry 
will  appear  in  the  command  queue  area. 

6.  Select  the  “Export"  tab.  This  will  cause  another  GUI  to  appear. 

7.  In  the  first  empty  field  type  the  name  of  the  exported  file. 

8.  Select  the  point  where  you  want  the  Fields  Calculator  to  compute  the  fields  in 
the  model.  There  are  two  ways  to  do  this.  First  is  to  have  HFSS  calculate  the 
locations  by  specifying  the  start  point,  end  points  and  spacing.  The  second 
way  is  to  create  a  tab  delimited  file  with  the  suffix  of  *.pts  and  import  the 
points.  The  file  should  have  no  column  header,  should  have  the  columns 
arranged  in  x,  y,  z  order  and  should  be  in  units  of  meters. 

9.  After  selecting  “ok”,  HFSS  will  export  a  file  with  the  field  data. 
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Near  Field  Line  Setup 


The  near  field  line  setup  is  an  option  in  HFSS  that  allows  the  user  to  solve  for  the  near 
fields  along  a  line.  It  can  be  any  polyline  of  one  or  more  segments.  It  is  necessary  to 
have  the  near  field  line  setup  to  measure  values  along  a  specific  line  such  as  a  radial 
line  leading  away  from  the  source.  Before  near  field  data  can  be  plotted  or  exported  a 
setup  must  be  defined.  This  is  done  in  the  “Radiation'’  option  in  the  project  manager. 

1 .  Select  a  line  or  lines  that  you  wish  to  use. 

2.  Select  “Radiation/Insert  Near  Field  Setup/Line.”  A  new  window  will  appear 
with  a  couple  options  to  select. 

3.  Specify  the  name  of  the  set-up. 

4.  Select  what  geometry  is  desired  for  the  set  up.  If  the  right  line  was  not  initially 
elected,  it  can  be  selected  now. 

5.  Input  the  number  of  evenly  spaced  points  that  will  be  calculated  along  the 
geometry. 

6.  All  other  options  should  be  left  as  default. 

This  process  can  also  be  used  to  get  points  to  use  in  the  fields  calculator.  By  right 
clicking  on  the  setup  for  the  geometry  of  interest,  select  “Compute  Max  Parameters.” 
A  new  window  appears  and  gives  several  options  to  choose,  in  particular  the  solution 
to  use  and  the  frequency.  Select  “ok”  to  continue.  HFSS  will  do  some  calculation  and 
another  window  will  appear  giving  a  summary  of  the  near  field  values  to  be  exported. 
Select  “Export  Fields”  and  a  *.csv  file  will  be  created  that  will  contain  points  along 
the  line  and  field  values.  This  file  can  then  be  opened  in  excel  and  edited  to  meet  the 
format  required  for  the  fields  calculator  by  deleting  the  headers,  index,  and  field 
values  so  that  all  that  is  left  is  the  x,  y,  z  points.  Save  the  file  as  a  tab  delimited  file  and 
change  the  file  suffix  to  “*.pts.” 

Validation 

The  following  discussion  hinges  on  the  validation  of  HFSS  data.  The  actual 
electromagnetic  phenomena  associated  w  ith  this  data  are  of  secondary  importance. 

•  Figures  1  and  2:  The  first  validation  scenario  considers  a  1,000  Hz  electric 
source  in  homogeneous,  lossy  water  with  a  conductivity  of  0.018  S/m.  Data 
generated  by  HFSS  and  by  the  standard  electric  dipole  solution  are  plotted  and 
compared.  The  correlation  between  data  sets  is  quite  good,  but  a  calibration 
multiplier  of  six  was  needed  in  the  HFSS  data.  The  reason  for  this  factor  is  not 
readily  apparent  to  us,  but  further  studies  performed  by  us  have  suggested  that 
the  factor  is  proportional  to  f  :.  Since  the  wave  number  is  also  proportional  to 
f1 2  in  lossy  media,  we  surmise  that  the  calibration  factor  is  related  to  wave 
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number. 

•  Figures  3  and  4:  The  second  validation  scenario  considers  a  100  Hz  source  in 
free  space.  Data  generated  by  HFSS  and  by  the  standard  electric  dipole 
solution  are  plotted  and  compared.  Although  the  electric  field  comparison  is 
quite  good,  the  magnetic  field  comparison  is  not.  The  trend  in  the  magnetic 
field  is  correct,  but  the  HFSS  data  is  sporadic  about  the  trend  line.  Since  HFSS 
is  an  electric  field  solver,  this  sporadic  data  is  associated  with  the 
post-processing  step  that  converts  electric  field  data  into  magnetic  field  data 
via  Ampere's  law.  The  quality  of  the  magnetic  field  data  is  improved  by  mesh 
refinement,  since  the  mesh  dictates  the  quality  of  the  numerical  derivatives.  In 
addition  to  this  issue,  a  multiplication  factor  of  sixteen  was  needed  to  get  the 
HFSS  data  to  agree  with  the  standard  dipole  solution.  The  need  for  this  factor 
is  still  unknown. 

•  The  data  of  Figures  6,  7,  8  and  9  are  associated  with  the  three 
layered-geometry  of  Figure  5.  The  conductivity  of  the  water  is  0.018  S/m;  the 
conductivity  of  the  floor  is  0.012  S/m.  For  Figures  6  and  7,  the  excitation  is  of 
the  electric  kind;  for  Figure  8  and  9,  the  excitation  is  of  the  magnetic  kind.  The 
operating  frequency  in  both  cases  is  1,000  Hz.  In  both  cases  the  source  and 
observation  points  are  located  1 1  m  below  the  surface  and  5  m  below  the 
surface,  respectively.  In  all  four  figures,  the  correlation  between  the  HFSS  data 
and  the  Sommerfeld  (WSU)  data  is  quite  good.  Again,  however,  a 
multiplication  factor  is  needed  in  the  HFSS  data  to  obtain  good  correlation.  In 
the  case  of  Figures  6  and  7,  the  factor  is  six;  in  the  case  of  Figures  8  and  9,  the 
factor  is  sixty.  Further  research  is  needed  to  understand  and  to  quantify  this 
factor. 
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HFSS  vs.  Analytical  Solution 


Homogenous  Water 


Figure  1 


HFSS  vs.  Analytical  Solution 

Homogenous  Water 


Figure  2 
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H-Field  Amplitude  (A/m) 


HFSS  vs.  Analytical  Solution 


Homogenous  Air 


Figure  3 


HFSS  vs.  Analytical  Solution 


Homogeneous  Air 


Figure  4 
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Figure  5:  A  depiction  of  a  flat,  three-layered  Lake. 
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HFSSvs.  WSU 

Three  Layer  Lake 
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Figure  7 


HFSS  vs.  Analytical  Solution 


Three  Layer  Lake 
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Figure  8 
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HFSSvs.  WSU 

Three  Layer  Lake 


Figure  9 


Future  Work 

It  is  clear  from  the  preceding  discussion  that  HFSS  is  indeed  a  viable  tool  for 
computing  ELF  signals  in  water.  However,  there  still  exists  several  issues  that  need  to 
be  resolved  to  make  this  tool  highly  effective.  These  are: 

•  Multiplication  Factor:  As  noted  in  the  previous  data,  a  frequency-dependent 
multiplication  factor  is  needed  to  get  correlation  between  the  HFSS  data  and 
the  closed-form  solutions.  The  reason  for  this  factor  is  not  understood  and 
needs  to  be  investigated. 

•  Electric  Field  Excitation,  Magnetic  Field  Observation:  When  an  ELF  electric 
source  is  used  to  excite  a  magnetic  field,  the  corresponding  magnetic  field  data 
in  HFSS  follows  the  correct  trend,  but  the  data  is  still  nevertheless  quite 
sporadic.  We  believe  that  this  issue  is  directly  related  to  the  mesh.  Some  kind 
of  mesh  refinement  is  believed  to  be  needed  to  resolve  this  problem. 

•  Meshing:  Given  that  the  Lake  has  regions  where  the  floor  and  surrounding 
landscape  spatially  varies  in  large  amounts  and  has  other  regions  where  such 
variations  are  little,  the  finite-element  mesh  needs  to  be  optimized  to  obtain 
both  rapid  and  accurate  solutions.  A  meshing  software  package  has  been 
purchased  to  assist  HFSS  with  the  meshing  process.  Use  of  this  package  will 
be  a  key  activity  of  Phase  Two. 
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•  Computer  Resources:  HFSS  is  a  computationally  intensive  program.  Its 
ability  to  solve  large,  complex  probelms  is  directly  linked  to  the  available 
computer  resources.  In  the  case  of  this  ELF  project,  a  resource  study  will  be 
conducted  to  ascertain  the  required  computational  resources  that  HFSS  will 
need  to  accomplish  the  stated  tasks. 
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Appendix  D:  Quasi-Static  Method 

Robert  Rebich  and  Christopher  Wagner 
Department  of  Electrical  and  Computer  Engineering,  University  of  Idaho 

Introduction 

For  ELF  frequencies,  the  corresponding  wavelengths  are  very  large  compared  to  the 
domain  size,  which  suggests  that  wavelike  motion  is  not  a  main  attribute  of  the 
electromagnetic  field.  Therefore  an  electrostatic  solution  can  be  applied  to  approximate 
the  fields.  For  a  system  to  be  regarded  as  quasi-static  the  time  rate  of  change  in  the 
electric  field  is  negligible  for  magnetic  source  excitation;  the  time  rate  of  change  in  the 
magnetic  field  is  negligible  for  electric  source  excitation.  This  allows  us  to  write 
electric  field  solutions  in  terms  of  a  scalar  potential  and  magnetic  field  solutions  in 
terms  of  a  vector  potential.  In  this  work,  we  have  only  concerned  ourselves  with 
quasi-electrostatic  (QES)  effects. 

The  QES  system  is  shown  in  Figure  I.  The  materials  that  define  the  two  layers  are 
characterized  by  their  conductivity  and  permittivity.  The  experiment  took  place  on  a 
lake;  therefore  the  two  material  media  model  will  be  air  and  water  with  their 
corresponding  material  characteristics  shown  in  the  figure.  The  green  dots  in  the  water 
represent  a  set  of  electrodes  with  opposite  charges  and  the  green  dots  in  the  air  represent 
the  images  electrodes.  The  image  electrodes  produce  the  same  fields  at  the  EMA  that 
the  reflected  fields  off  the  interface  from  the  original  electrodes.  The  electrodes  are 
located  at  a  vertical  distance  h  from  the  interface  and  separated  by  a  distance  d.  The  red 
dot  represents  the  EMA  (x,  y,  z).  The  model  replicates  the  experiment  by  moving  the 
electrodes  along  a  radial  axis  and  measuring  the  fields  at  the  EMA. 

To  energize  the  electrodes,  a  voltage  V  is  impressed  between  them.  The  experimental 
voltage  between  the  electrodes  is  unknown,  but  the  current  was  recorded.  To  find  the 
voltage,  we  estimated  the  impedance  between  the  electrodes  and  multiplied  that 
impedance  by  the  current. 
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Figure  1:  Geometrical  rendition  of  the  quasi-static  geometry. 

The  voltage  can  be  determined  anywhere  in  the  domain  by  a  simple  electrostatics 
relationship.  The  total  voltage  is  the  superposition  of  each  contribution  made  by  each  of 
the  real  and  image  electrodes.  The  voltages  in  the  water  and  air  are  given  as, 

V  =  — ^ L  Z>0 

iK(at+jcjs1)  Lftt  rtI  r2t  r.2J 


v  =  Jzi _ _  [21  _  JL 

4k  (ert+/cj  jrt)  Lftl  r21J ' 


Z<() 


Here  R  is  the  reflection  coefficient  at  the  interface  and  Tis  the  transmission  coefficient 
at  the  interface.  The  posistion  vectors  are  functions  of  x,  y  and  z  and  are  as  follows: 


rn 


ra  =  JpI  +  Cz  +  -l): 
r«  =  Jpf+b  +  zJ* 


pL=  VCx- xL):  +  (y  -  yL}2 

Pi  =  Vu  -x.)-  +  (y  -y2> 


The  subscripted  coordinates  refer  to  either  the  first  charge  (left)  or  the  second  charge 
(right).  Once  the  position  vectors  and  potentials  are  known,  the  total  electric  field  can 
be  determined  through  the  gradient  of  the  voltage: 

E  =  —  7V 

The  validity  of  the  code  is  based  upon  the  assumption  that  the  charge  on  the  spheres  is 
uniformly  distributed.  This  assumption  will  fail  when  the  two  spheres  are  close  to  each 
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other  or  when  the  spheres  are  close  to  the  material  interface. 


A  comparison  of  the  QES  solution  versus  experimental  and  other  solution  methods  is 
shown  in  Figure  2.  The  dipole  for  the  experiment  was  oriented  down  the  positive  x-axis, 
which  resulted  in  a  very  weak  component  of  Ev\  therefore  only  Ex  (solid  curves)  and  E-_ 
(dashed  curves)  are  represented.  The  plot  shows  very  good  comparisons  for  the  QES 
and  Sommerfeld  (WSU)  solutions.  After  about  100  m  the  WSU  Ez  component  flares 
away  from  the  QES  E:.  This  is  expected  since  the  QES  solution  cannot  account  for 
wave-like  effects,  unlike  the  WSU  solution.  The  experimental  data  from  run  #1 1 17  are 
also  comparable  to  known  analytical  solutions.  The  characteristics  of  the  curves 
matched  the  characteristics  of  the  analytical  solutions. 


Scenario  #5  -  Electric  Hull  - 1 00  Hz  -  0.5  A*m 

E x  =  solid  E  z  =  dashed  (Run  1117) 


Figure  2:  Field  comparison  for  Run  #1117.  Source  is  0.5  A-m. 
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